#Table A10

#directory
#set working directory to the path PvP_Replication
#setwd("~/PvP_Replication")

#packages
library(tidyverse) #for data cleaning and manipulation
library(estimatr) #estimation
library(texreg) #latex tables

#load main dataset
pvp_data <- read.csv("PvP_data/PvP_data_main.csv")

#transform independent variable
#create log + 1 cultivation variable
pvp_data <- pvp_data %>% 
  mutate(logcult = log(maxcult + 1))

#subset data to municipalities with prior FARC presence
pvp_farc_subset <- pvp_data %>% filter(farc_presence == 1)

#prep weather and soil variables for iv
pvp_farc_subset <- pvp_farc_subset %>% mutate(
  draindummy = as.factor(ifelse(DRAIN == "M", 1, 0)), #moderate drainage
  humiddist = case_when(humid_mean < 85 ~ (85 - humid_mean)^2, humid_mean > 85 ~ (humid_mean - 85)^2, TRUE ~ 0 ), #distance from optimal humidity point
  temprangedist = case_when(maxtemp > 27 ~ (maxtemp - 27)^2, mintemp < 14  ~ (14 - mintemp)^2, TRUE ~ 0), #distance from optimal temp, wide band
)

#instrument for continuous independent variable
pvp_farc_subset$ols_inst <- predict(lm(logcult ~ PHAQ + temprangedist + TOTN + TOTC + sun_mean + humiddist + raintot_me + draindummy, pvp_farc_subset))

#scale variables to facilitate comparisons

pvp_farc_subset <- pvp_farc_subset %>% 
  mutate(rugged_scale = as.numeric(scale(rugged_mean)), 
                                   log_cult_scaled = as.numeric(scale(logcult)))

#estimate models
ols_altrug_only <- lm_robust(dissident_presence ~ rugged_scale, pvp_farc_subset)
ols_altrug_full <- lm_robust(dissident_presence ~ log_cult_scaled + rugged_scale, pvp_farc_subset)
iv_altrug_full <- iv_robust(dissident_presence ~ log_cult_scaled + rugged_scale | ols_inst + rugged_scale, data = pvp_farc_subset)

#generate table 
texreg(
  list(ols_altrug_only, iv_altrug_full, ols_altrug_full),
  custom.coef.names	= c(NA, "Terrain Ruggedness (TRI)", "Coca Cultivation"),
  digits = 2, include.ci = FALSE, single.row = FALSE, include.fstatistic = TRUE, 
  include.rmse = FALSE, include.rsquared = FALSE, include.adjrs = FALSE, 
  include.nobs = TRUE, stars = c(0.001, 0.01, 0.05), float.pos = "h", 
  custom.header = list("OLS" = 1, "IV" = 2, "OLS" = 3),
  caption.above	= TRUE, caption = "Alternative Measure of Terrain Ruggedness")